homework 4, version 1
Submission by: Jazzy Doe (jazz@mit.edu)
Homework 4: Epidemic modeling I
18.S191, fall 2020
This notebook contains built-in, live answer checks! In some exercises you will see a coloured box, which runs a test case on your code, and provides feedback based on the result. Simply edit the code, run it, and the check runs again.
For MIT students: there will also be some additional (secret) test cases that will be run as part of the grading process, and we will look at your notebook and write comments.
Feel free to ask questions!
"Jazzy Doe"
"jazz"
xxxxxxxxxx# edit the code below to set your name and kerberos ID (i.e. email without @mit.edu)student = (name = "Jazzy Doe", kerberos_id = "jazz")# you might need to wait until all other cells in this notebook have completed running. # scroll around the page to see what's upLet's create a package environment:
xxxxxxxxxxbegin using Pkg Pkg.activate(mktempdir())endxxxxxxxxxxbegin Pkg.add(["Plots", "PlutoUI",]) using Plots, Statistics plotly() using PlutoUIendExercise 1: Modelling recovery
In this exercise we will investigate a simple stochastic (probabilistic) model of recovery from an infection and the time
In this model, an individual who is infected has a constant probability
Exercise 1.1 - Probability distributions
👉 Define the function bernoulli(p), which returns true with probability false with probability
bernoulli (generic function with 1 method)xxxxxxxxxxfunction bernoulli(p::Number) return rand() <= pendGot it!
Keep it up!
👉 Write a function recovery_time(p) that returns the time taken until the person recovers.
recovery_time (generic function with 1 method)x
function recovery_time(p) if p ≤ 0 throw(ArgumentError("p must be positive: p = 0 cannot result in a recovery")) end recovered = false days = 1 while recovered == false if bernoulli(p) recovered = true else days += 1 end end # Your code here. See the comment below about the p ≤ 0 case. return daysendGot it!
Great!
Hint
Remember to always re-use work you have done previously: in this case you should re-use the function bernoulli.
We should always be aware of special cases (sometimes called "boundary conditions"). Make sure not to run the code with ArgumentError as follows:
throw(ArgumentError("..."))
with a suitable error message.
👉 What happens for
immediate recovery
xxxxxxxxxxinterpretation_of_p_equals_one = md"""immediate recovery """Exercise 1.2
👉 Write a function do_experiment(p, N) that runs the function recovery_time N times and collects the results into a vector.
do_experiment (generic function with 1 method)xxxxxxxxxxfunction do_experiment(p, N) results = zeros(N) for i ∈ 1:N results[i] = recovery_time(p) end return resultsend3.0
1.0
1.0
1.0
2.0
1.0
2.0
5.0
2.0
1.0
1.0
1.0
3.0
2.0
2.0
2.0
3.0
1.0
2.0
1.0
xxxxxxxxxxsmall_experiment = do_experiment(0.5, 20)Exercise 1.3
👉 Write a function frequencies(data) that calculates and returns the frequencies (i.e. probability distribution) of input data.
The input will be an array of integers, with duplicates, and the result will be a dictionary that maps each occured value to its frequency in the data.
For example,
frequencies([7, 8, 9, 7])
should give
Dict(
7 => 0.5,
8 => 0.25,
9 => 0.25
)
As with any probability distribution, it should be normalised to
frequencies (generic function with 1 method)xxxxxxxxxxbeginusing StatsBase function frequencies(values) return StatsBase.countmap(values) endend2.0
7
3.0
3
5.0
1
1.0
9
xxxxxxxxxxfrequencies(small_experiment)Hint
Do you remember how we worked with dictionaries in Homework 3? You can create an empty dictionary using Dict(). You may want to use either the function haskey or the function get on your dictionary – check the documentation for how to use these functions.
Let's run an experiment with
1.0
5.0
3.0
2.0
1.0
6.0
4.0
7.0
7.0
6.0
2.0
11.0
5.0
2.0
2.0
4.0
1.0
1.0
4.0
3.0
4.0
1.0
3.0
2.0
2.0
2.0
6.0
2.0
2.0
4.0
3.0
17.0
10.0
7.0
3.0
12.0
7.0
5.0
3.0
3.0
1.0
1.0
1.0
6.0
3.0
6.0
5.0
1.0
2.0
10.0
xxxxxxxxxxlarge_experiment = do_experiment(0.25, 10_000) # (10_000 is just 10000 but easier to read)The frequencies dictionary is difficult to interpret on its own, so instead, we will plot it, i.e. plot
Plots.jl comes with a function bar, which does exactly what we want:
xxxxxxxxxxbar(frequencies(large_experiment))Great! Feel free to experiment with this function, try giving it a different array as argument. Plots.jl is pretty clever, it even works with an array of strings!
Exercise 1.4
Next, we want to add a new element to our plot: a vertical line. To demonstrate how this works, here we added a vertical line at the maximum value.
To write this function, we first create a base plot, we then modify that plot to add the vertical line, and finally, we return the plot. More on this in the next info box.
frequencies_plot_with_maximum (generic function with 1 method)xxxxxxxxxxfunction frequencies_plot_with_maximum(data::Vector) base = bar(frequencies(data)) vline!(base, [maximum(data)], label="maximum") return baseendxxxxxxxxxxfrequencies_plot_with_maximum(large_experiment)Note about plotting
Plots.jl has an interesting property: a plot is an object, not an action. Functions like
plot,bar,histogramdon't draw anything on your screen - they just return aPlots.Plot. This is a struct that contains the description of a plot (what data should be plotted in what way?), not the picture.So a Pluto cell with a single line,
plot(1:10), will show a plot, because the result of the functionplotis aPlotobject, and Pluto just shows the result of a cell.Modifying plots
Nice plots are often formed by overlaying multiple plots. In Plots.jl, this is done using the modifying functions:
plot!,bar!,vline!, etc. These take an extra (first) argument: a previous plot to modify.For example, to plot the
sin,cosandtanfunctions in the same view, we do:function sin_cos_plot() T = -1.0:0.01:1.0 result = plot(T, sin.(T)) plot!(result, T, cos.(T)) plot!(result, T, tan.(T)) return result end💡 This example demonstrates a useful pattern to combine plots:
Create a new plot and store it in a variable
Modify that plot to add more elements
Return the plot
It is recommended that these 3 steps happen within a single cell. This can prevent some strange glitches when re-running cells. There are three ways to group expressions together into a single cell:
begin,letandfunction. More on this later!
👉 Write the function frequencies_plot_with_mean that calculates the mean recovery time and displays it using a vertical line.
frequencies_plot_with_mean (generic function with 1 method)x
function frequencies_plot_with_mean(data) # start out by copying the frequencies_plot_with_maximum function base = bar(frequencies(data)) vline!(base, [mean(data)], label="mean") return baseendxxxxxxxxxxfrequencies_plot_with_mean(large_experiment)👉 Write an interactive visualization that draws the histogram and mean for p_interactive and N_interactive, instead of just p and N.
xxxxxxxxxx p_interactive Slider(0.01:.01:1; default=0.25, show_value=true)xxxxxxxxxx N_interactive Slider(1:10_000, default=5000, show_value=true)xxxxxxxxxxfrequencies_plot_with_mean(do_experiment(p_interactive, N_interactive))As you separately vary
xxxxxxxxxxletλ = 2.5p = plot()x1 = 1:.01:10; y1 = (λ*ℯ).^-(λ*x1);plot!(p, x1, y1)endExercise 1.5
👉 What shape does the distribution seem to have? Can you verify that by adding a second plot with the expected shape?
Closely resembles exponential distribution
x
md"""Closely resembles exponential distribution"""Exercise 1.6
👉 Use
begin n = 10_000 xs =0.001:0.01:1 results = [mean(do_experiment(p, n)) for p ∈ xs] plot(xs, results, title="Avg Days to Recovery") ylabel!("Avg Days to Recover") xlabel!("Probability of Recovery")endExercise 2: Agent-based model for an epidemic outbreak – types
In this and the following exercises we will develop a simple stochastic model for combined infection and recovery in a population, which may exhibit an epidemic outbreak (i.e. a large spike in the number of infectious people). The population is well mixed, i.e. everyone is in contact with everyone else. [An example of this would be a small school or university in which people are constantly moving around and interacting with each other.]
The model is an individual-based or agent-based model: we explicitly keep track of each individual, or agent, in the population and their infection status. For the moment we will not keep track of their position in space; we will just assume that there is some mechanism, not included in the model, by which they interact with other individuals.
Exercise 2.1
Each agent will have its own internal state, modelling its infection status, namely "susceptible", "infectious" or "recovered". We would like to code these as values S, I and R, respectively. One way to do this is using an enumerated type or enum. Variables of this type can take only a pre-defined set of values; the Julia syntax is as follows:
xxxxxxxxxx InfectionStatus S I RWe have just defined a new type InfectionStatus, as well as names S, I and R that are the (only) possible values that a variable of this type can take.
👉 Define a variable test_status whose value is S.
S::InfectionStatus = 0xxxxxxxxxxtest_status = S👉 Use the typeof function to find the type of test_status.
Enum InfectionStatus:
S = 0
I = 1
R = 2xxxxxxxxxxtypeof(test_status)👉 Convert x to an integer using the Integer function. What value does it have? What values do I and R have?
0
1
2
xxxxxxxxxxInt(S), Int(I), Int(R)Exercise 2.2
For each agent we want to keep track of its infection status and the number of other agents that it infects during the simulation. A good solution for this is to define a new type Agent to hold all of the information for one agent, as follows:
Agentxxxxxxxxxxbegin mutable struct Agent status::InfectionStatus num_infected::Int64 end Agent() = Agent(S,0)endWhen you define a new type like this, Julia automatically defines one or more constructors, which are methods of a generic function with the same name as the type. These are used to create objects of that type.
👉 Use the methods function to check how many constructors are pre-defined for the Agent type.
- Agent() in Main.workspace3 at /Users/Rich/CSprojects/Tutorials/18.S191CompThinking/hw4.jl#==#ae4ac4b4-041f-11eb-14f5-1bcde35d18f2:7
- Agent(status::Main.workspace3.InfectionStatus, num_infected::Int64) in Main.workspace3 at /Users/Rich/CSprojects/Tutorials/18.S191CompThinking/hw4.jl#==#ae4ac4b4-041f-11eb-14f5-1bcde35d18f2:3
- Agent(status, num_infected) in Main.workspace3 at /Users/Rich/CSprojects/Tutorials/18.S191CompThinking/hw4.jl#==#ae4ac4b4-041f-11eb-14f5-1bcde35d18f2:3
xxxxxxxxxxmethods(Agent)👉 Create an agent test_agent with status S and num_infected equal to 0.
S::InfectionStatus = 0
0
xxxxxxxxxxtest_agent = Agent(S, 0)👉 For convenience, define a new constructor (i.e. a new method for the function) that takes no arguments and creates an Agent with status S and number infected 0, by calling one of the default constructors that Julia creates. This new method lives outside (not inside) the definition of the struct. (It is called an outer constructor.)
(In Pluto, multiple methods for the same function need to be combined in a single cell using a begin end block.)
Let's check that the new method works correctly. How many methods does the constructor have now?
S::InfectionStatus = 0
0
xxxxxxxxxxAgent()Exercise 2.3
👉 Write functions set_status!(a) and set_num_infected!(a) which modify the respective fields of an Agent. Check that they work. [Note the bang ("!") at the end of the function names to signify that these functions modify their argument.]
set_num_infected! (generic function with 1 method)xxxxxxxxxxbegin function set_status!(agent::Agent, new_status::InfectionStatus) agent.status = new_status end function set_num_infected!(agent::Agent, num_infected::Int64) agent.num_infected = num_infected endendGot it!
Great!
👉 We will also need functions is_susceptible and is_infected that check if a given agent is in those respective states.
is_susceptible (generic function with 1 method)xxxxxxxxxxfunction is_susceptible(agent::Agent) return agent.status == Sendis_infected (generic function with 1 method)xxxxxxxxxxfunction is_infected(agent::Agent) return agent.status == IendGot it!
Keep it up!
Exericse 2.4
👉 Write a function generate_agents(N) that returns a vector of N freshly created Agents. They should all be initially susceptible, except one, chosen at random (i.e. uniformly), who is infectious.
Main.workspace3.generate_agentsxxxxxxxxxx"""Generates N agents with 1 random agent infected"""function generate_agents(N::Integer) infected_idx = rand(1:N) agents = [Agent() for i in 1:N] agents[infected_idx].status = I return agentsendGot it!
Fantastic!
We will also need types representing different infections.
Let's define an (immutable) struct called InfectionRecovery with parameters p_infection and p_recovery. We will make it a subtype of an abstract AbstractInfection type, because we will define more infection types later.
xxxxxxxxxxabstract type AbstractInfection endxxxxxxxxxxmutable struct InfectionRecovery <: AbstractInfection p_infection p_recoveryendExercise 2.5
👉 Write a function interact! that takes an affected agent of type Agent, an source of type Agent and an infection of type InfectionRecovery. It implements a single (one-sided) interaction between two agents:
If the
agentis susceptible and thesourceis infectious, then thesourceinfects ouragentwith the given infection probability. If thesourcesuccessfully infects the other agent, then itsnum_infectedrecord must be updated.If the
agentis infected then it recovers with the relevant probability.Otherwise, nothing happens.
interact! (generic function with 2 methods)xxxxxxxxxxbegin function interact!(agent::Agent, source::Agent, infection::InfectionRecovery) if is_susceptible(agent) && is_infected(source) if bernoulli(infection.p_infection) set_status!(agent, I) set_num_infected!(source, source.num_infected + 1) end elseif is_infected(agent) if bernoulli(infection.p_recovery) set_status!(agent, R) end end end function interact!(agent::Agent, source::Agent, infection::Reinfection) if is_susceptible(agent) && is_infected(source) if bernoulli(infection.p_infection) set_status!(agent, I) set_num_infected!(source, source.num_infected + 1) end elseif is_infected(agent) if bernoulli(infection.p_recovery) set_status!(agent, R) end elseif agent.status == R set_status!(agent, S) end endendPlay around with the test case below to test your function! Try changing the definitions of agent, source and infection. Since we are working with randomness, you might want to run the cell multiple times.
xxxxxxxxxxmd"""Play around with the test case below to test your function! Try changing the definitions of `agent`, `source` and `infection`. Since we are working with randomness, you might want to run the cell multiple times."""I::InfectionStatus = 1
0
I::InfectionStatus = 1
1
xxxxxxxxxxlet agent = Agent(S, 0) source = Agent(I, 0) infection = InfectionRecovery(0.9, 0.5) interact!(agent, source, infection) (agent=agent, source=source)endGot it!
Your function treats the susceptible agent case correctly!
Got it!
Your function treats the infectious agent case correctly!
Got it!
Your function treats the recovered agent case correctly!
Exercise 3: Agent-based model for an epidemic outbreak – Monte Carlo simulation
In this exercise we will build on Exercise 2 to write a Monte Carlo simulation of how an infection propagates in a population.
Make sure to re-use the functions that we have already written, and introduce new ones if they are helpful! Short functions make it easier to understand what the function does and build up new functionality piece by piece.
You should not use any global variables inside the functions: Each function must accept as arguments all the information it requires to carry out its task. You need to think carefully about what the information each function requires.
Exercise 3.1
👉 Write a function step! that takes a vector of Agents and an infection of type InfectionRecovery. It implements a single step of the infection dynamics as follows:
Choose two random agents: an
agentand asource.Apply
interact!(agent, source, infection).Return
agents.
step! (generic function with 1 method)x
function step!(agents::Vector{Agent}, infection::AbstractInfection) n_agents = length(agents) a1_idx, a2_idx = 1, 1 while a1_idx == a2_idx # need to insure agent1 != agent2 a1_idx, a2_idx = rand(1:n_agents), rand(1:n_agents) end a1, a2 = agents[a1_idx], agents[a2_idx] interact!(a1, a2, infection) return agentsend👉 Write a function sweep!. It runs step!
sweep! (generic function with 1 method)function sweep!(agents::Vector{Agent}, infection::AbstractInfection) n_agents = length(agents) for i in 1:n_agents step!(agents, infection) endend👉 Write a function simulation that does the following:
Generate the
agents.Run
sweep!a number of times. Calculate and store the total number of agents with each status at each step in variablesS_counts,I_countsandR_counts.Return the vectors
S_counts,I_countsandR_countsin a named tuple, with keysS,IandR.
You've seen an example of named tuples before: the student variable at the top of the notebook!
Feel free to store the counts in a different way, as long as the return type is the same.
simulation (generic function with 1 method)function simulation(N::Integer, T::Integer, infection::AbstractInfection) agents = generate_agents(N) S_counts, I_counts, R_counts = zeros(T), zeros(T), zeros(T) for t in 1:T sweep!(agents, infection) for a in agents if is_susceptible(a) S_counts[t] += 1 elseif is_infected(a) I_counts[t] += 1 else R_counts[t] += 1 end end end return (S=S_counts, I=I_counts, R=R_counts)end2.0
2.0
2.0
2.0
2.0
2.0
2.0
2.0
2.0
2.0
2.0
2.0
2.0
2.0
2.0
2.0
2.0
2.0
2.0
2.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
xxxxxxxxxxsimulation(3, 20, InfectionRecovery(0.9, 0.2))99.0
99.0
99.0
99.0
99.0
99.0
99.0
99.0
99.0
99.0
99.0
99.0
99.0
99.0
99.0
99.0
99.0
99.0
99.0
99.0
99.0
99.0
99.0
99.0
99.0
99.0
99.0
99.0
99.0
99.0
99.0
99.0
99.0
99.0
99.0
99.0
99.0
99.0
99.0
99.0
99.0
99.0
99.0
99.0
99.0
99.0
99.0
99.0
99.0
99.0
1.0
1.0
1.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
xxxxxxxxxxsimulation(100, 1000, InfectionRecovery(0.005, 0.2))let run_basic_sir N = 100 T = 1000 sim = simulation(N, T, InfectionRecovery(0.02, 0.002)) result = plot(1:T, sim.S, ylim=(0, N), label="Susceptible") plot!(result, 1:T, sim.I, ylim=(0, N), label="Infectious") plot!(result, 1:T, sim.R, ylim=(0, N), label="Recovered")endWe used a let block in this cell to group multiple expressions together, but how is it different from begin or function?
function vs. begin vs. let
Writing functions is a way to group multiple expressions (i.e. lines of code) together into a mini-program. Note the following about functions:
A function always returns one object.[1] This object can be given explicitly by writing
return x, or implicitly: Julia functions always return the result of the last expression by default. Sof(x) = x+2is the same asf(x) = return x+2.Variables defined inside a function are not accessible outside the function. We say that function bodies have a local scope. This helps to keep your program easy to read and write: if you define a local variable, then you don't need to worry about it in the rest of the notebook.
There are two other ways to group epxressions together that you might have seen before:
beginandlet.begin
beginwill group expressions together, and it takes the value of its last subexpression.We use it in this notebook when we want multiple expressions to always run together.
let
letalso groups multiple expressions together into one, but variables defined inside of it are local: they don't affect code outside of the block. So likebegin, it is just a block of code, but likefunction, it has a local variable scope.We use it when we want to define some local (temporary) variables to produce a complicated result, without interfering with other cells. Pluto allows only one definition per global variable of the same name, but you can define local variables with the same names whenever you wish!
1
Even a function like
f(x) = returnreturns one object: the object
nothing— try it out!
Exercise 3.2
Alright! Every time that we run the simulation, we get slightly different results, because it is based on randomness. By running the simulation a number of times, you start to get an idea of the mean behaviour of our model. This is the essence of a Monte Carlo method! You use computer-generated randomness to generate samples.
Instead of pressing the button many times, let's have the computer repeat the simulation. In the next cells, we run your simulation num_simulations=20 times with
Every single simulation returns a named tuple with the status counts, so the result of multiple simulations will be an array of those. Have a look inside the result, simulations, and make sure that its structure is clear.
repeat_simulations (generic function with 1 method)x
function repeat_simulations(N, T, infection, num_simulations) N = 100 T = 1000 map(1:num_simulations) do _ simulation(N, T, infection) endend99.0
99.0
99.0
99.0
99.0
99.0
99.0
99.0
99.0
99.0
99.0
99.0
99.0
99.0
99.0
98.0
97.0
97.0
97.0
97.0
97.0
97.0
97.0
97.0
97.0
97.0
96.0
96.0
96.0
96.0
96.0
96.0
96.0
96.0
96.0
96.0
96.0
95.0
95.0
95.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
2.0
3.0
3.0
3.0
3.0
3.0
3.0
3.0
3.0
3.0
3.0
4.0
4.0
4.0
4.0
4.0
4.0
4.0
4.0
4.0
4.0
4.0
5.0
5.0
5.0
28.0
28.0
28.0
28.0
27.0
27.0
27.0
27.0
27.0
27.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
72.0
72.0
72.0
72.0
73.0
73.0
73.0
73.0
73.0
73.0
99.0
99.0
99.0
99.0
99.0
99.0
99.0
99.0
99.0
99.0
99.0
99.0
99.0
99.0
99.0
99.0
99.0
99.0
99.0
99.0
99.0
99.0
99.0
99.0
99.0
99.0
99.0
99.0
99.0
99.0
99.0
99.0
99.0
99.0
99.0
99.0
99.0
99.0
99.0
99.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
23.0
23.0
23.0
23.0
23.0
23.0
23.0
23.0
23.0
23.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
77.0
77.0
77.0
77.0
77.0
77.0
77.0
77.0
77.0
77.0
99.0
99.0
99.0
99.0
99.0
99.0
99.0
99.0
99.0
99.0
99.0
99.0
99.0
99.0
99.0
99.0
99.0
99.0
99.0
99.0
99.0
99.0
99.0
99.0
99.0
99.0
99.0
99.0
99.0
99.0
99.0
99.0
99.0
99.0
99.0
99.0
99.0
99.0
99.0
99.0
99.0
99.0
99.0
99.0
99.0
99.0
99.0
99.0
99.0
99.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
99.0
99.0
99.0
99.0
99.0
99.0
99.0
99.0
99.0
99.0
99.0
99.0
99.0
99.0
99.0
99.0
99.0
98.0
98.0
98.0
98.0
98.0
98.0
98.0
98.0
98.0
97.0
97.0
97.0
97.0
97.0
97.0
97.0
97.0
97.0
97.0
97.0
97.0
97.0
97.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
2.0
2.0
2.0
2.0
2.0
2.0
2.0
2.0
2.0
3.0
3.0
3.0
3.0
3.0
3.0
3.0
3.0
3.0
3.0
3.0
3.0
3.0
3.0
30.0
30.0
30.0
30.0
30.0
29.0
29.0
29.0
29.0
29.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
70.0
70.0
70.0
70.0
70.0
71.0
71.0
71.0
71.0
71.0
99.0
99.0
99.0
98.0
98.0
98.0
98.0
98.0
98.0
98.0
98.0
98.0
98.0
98.0
98.0
98.0
98.0
98.0
98.0
97.0
97.0
97.0
97.0
97.0
97.0
97.0
97.0
97.0
97.0
97.0
97.0
97.0
96.0
96.0
96.0
96.0
96.0
96.0
96.0
96.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
1.0
1.0
1.0
2.0
2.0
2.0
2.0
2.0
2.0
2.0
2.0
2.0
2.0
2.0
2.0
2.0
2.0
2.0
2.0
3.0
3.0
3.0
3.0
3.0
3.0
3.0
3.0
3.0
3.0
3.0
3.0
3.0
4.0
4.0
4.0
4.0
4.0
4.0
4.0
4.0
18.0
18.0
18.0
18.0
17.0
17.0
17.0
17.0
17.0
17.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
82.0
82.0
82.0
82.0
83.0
83.0
83.0
83.0
83.0
83.0
99.0
99.0
99.0
99.0
99.0
99.0
99.0
99.0
99.0
99.0
99.0
99.0
99.0
99.0
99.0
99.0
99.0
99.0
99.0
99.0
99.0
99.0
99.0
99.0
99.0
99.0
99.0
99.0
99.0
99.0
99.0
99.0
99.0
99.0
99.0
99.0
99.0
99.0
99.0
99.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
28.0
28.0
28.0
28.0
28.0
28.0
28.0
28.0
28.0
28.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
72.0
72.0
72.0
72.0
72.0
72.0
72.0
72.0
72.0
72.0
99.0
98.0
98.0
98.0
98.0
98.0
98.0
98.0
98.0
98.0
98.0
98.0
98.0
98.0
98.0
98.0
98.0
98.0
98.0
98.0
98.0
98.0
98.0
98.0
98.0
98.0
98.0
97.0
97.0
97.0
97.0
96.0
96.0
96.0
96.0
96.0
96.0
96.0
96.0
95.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
1.0
2.0
2.0
2.0
2.0
2.0
2.0
2.0
2.0
2.0
2.0
2.0
2.0
2.0
2.0
2.0
2.0
2.0
2.0
2.0
2.0
2.0
2.0
2.0
2.0
2.0
2.0
3.0
3.0
3.0
3.0
4.0
4.0
4.0
4.0
4.0
4.0
4.0
4.0
5.0
23.0
23.0
23.0
23.0
23.0
23.0
23.0
23.0
23.0
23.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
77.0
77.0
77.0
77.0
77.0
77.0
77.0
77.0
77.0
77.0
99.0
99.0
99.0
99.0
99.0
99.0
99.0
99.0
99.0
99.0
99.0
99.0
99.0
99.0
99.0
99.0
99.0
99.0
99.0
99.0
99.0
99.0
99.0
99.0
99.0
99.0
99.0
99.0
99.0
99.0
99.0
99.0
99.0
99.0
99.0
99.0
99.0
99.0
99.0
99.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
31.0
31.0
31.0
31.0
31.0
31.0
31.0
31.0
31.0
31.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
69.0
69.0
69.0
69.0
69.0
69.0
69.0
69.0
69.0
69.0
99.0
99.0
99.0
99.0
99.0
99.0
99.0
99.0
99.0
99.0
99.0
99.0
98.0
98.0
98.0
98.0
98.0
98.0
98.0
97.0
97.0
97.0
97.0
97.0
97.0
97.0
97.0
97.0
97.0
97.0
97.0
97.0
97.0
97.0
97.0
97.0
97.0
97.0
97.0
97.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
2.0
2.0
2.0
2.0
2.0
2.0
2.0
3.0
3.0
3.0
3.0
3.0
3.0
3.0
3.0
3.0
3.0
3.0
3.0
3.0
3.0
3.0
3.0
3.0
3.0
3.0
3.0
3.0
21.0
21.0
21.0
21.0
21.0
21.0
21.0
21.0
21.0
21.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
78.0
78.0
78.0
78.0
78.0
78.0
78.0
78.0
78.0
78.0
99.0
99.0
99.0
99.0
99.0
99.0
99.0
99.0
99.0
99.0
99.0
98.0
98.0
98.0
98.0
98.0
98.0
98.0
98.0
98.0
97.0
97.0
97.0
97.0
97.0
97.0
97.0
97.0
97.0
97.0
97.0
97.0
97.0
97.0
97.0
97.0
97.0
97.0
97.0
97.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
2.0
2.0
2.0
2.0
2.0
2.0
2.0
2.0
2.0
3.0
3.0
3.0
3.0
3.0
3.0
3.0
3.0
3.0
3.0
3.0
3.0
3.0
3.0
3.0
3.0
3.0
3.0
3.0
3.0
21.0
20.0
20.0
20.0
20.0
20.0
20.0
20.0
20.0
20.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
79.0
80.0
80.0
80.0
80.0
80.0
80.0
80.0
80.0
80.0
99.0
99.0
98.0
98.0
98.0
98.0
98.0
97.0
97.0
97.0
97.0
97.0
96.0
96.0
96.0
96.0
96.0
96.0
96.0
96.0
96.0
96.0
96.0
96.0
96.0
96.0
96.0
95.0
95.0
95.0
95.0
95.0
95.0
95.0
95.0
95.0
95.0
95.0
95.0
95.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
1.0
1.0
2.0
2.0
2.0
2.0
2.0
3.0
3.0
3.0
3.0
3.0
4.0
4.0
4.0
4.0
4.0
4.0
4.0
4.0
4.0
4.0
4.0
4.0
4.0
4.0
4.0
5.0
5.0
5.0
5.0
5.0
5.0
5.0
5.0
5.0
5.0
5.0
5.0
5.0
14.0
14.0
14.0
14.0
14.0
14.0
14.0
14.0
14.0
14.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
86.0
86.0
86.0
86.0
86.0
86.0
86.0
86.0
86.0
86.0
99.0
99.0
99.0
99.0
99.0
99.0
99.0
99.0
99.0
99.0
99.0
99.0
99.0
99.0
99.0
99.0
99.0
99.0
99.0
99.0
99.0
99.0
99.0
99.0
99.0
99.0
99.0
99.0
99.0
99.0
99.0
99.0
99.0
99.0
99.0
99.0
99.0
99.0
99.0
99.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
26.0
26.0
26.0
26.0
26.0
25.0
25.0
25.0
25.0
24.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
74.0
74.0
74.0
74.0
74.0
75.0
75.0
75.0
75.0
76.0
99.0
99.0
99.0
99.0
99.0
99.0
99.0
99.0
99.0
99.0
99.0
99.0
99.0
99.0
99.0
99.0
99.0
99.0
99.0
99.0
98.0
98.0
98.0
98.0
98.0
98.0
98.0
98.0
98.0
98.0
98.0
98.0
98.0
98.0
98.0
98.0
98.0
98.0
98.0
98.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
2.0
2.0
2.0
2.0
2.0
2.0
2.0
2.0
2.0
2.0
2.0
2.0
2.0
2.0
2.0
2.0
2.0
2.0
2.0
2.0
29.0
29.0
29.0
28.0
28.0
28.0
28.0
28.0
28.0
28.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
71.0
71.0
71.0
72.0
72.0
72.0
72.0
72.0
72.0
72.0
99.0
99.0
99.0
99.0
99.0
99.0
99.0
99.0
99.0
99.0
99.0
99.0
99.0
99.0
99.0
99.0
99.0
99.0
99.0
99.0
99.0
99.0
99.0
99.0
99.0
99.0
99.0
99.0
99.0
99.0
99.0
99.0
99.0
99.0
99.0
99.0
99.0
99.0
99.0
99.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
33.0
32.0
32.0
32.0
32.0
32.0
32.0
32.0
32.0
32.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
67.0
68.0
68.0
68.0
68.0
68.0
68.0
68.0
68.0
68.0
99.0
99.0
99.0
99.0
99.0
99.0
99.0
99.0
99.0
99.0
99.0
99.0
99.0
99.0
99.0
99.0
99.0
99.0
99.0
99.0
98.0
98.0
98.0
98.0
98.0
98.0
98.0
98.0
98.0
98.0
98.0
98.0
98.0
98.0
98.0
98.0
98.0
98.0
98.0
98.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
2.0
2.0
2.0
2.0
2.0
2.0
2.0
2.0
2.0
2.0
2.0
2.0
2.0
2.0
2.0
2.0
2.0
2.0
2.0
2.0
26.0
25.0
25.0
25.0
25.0
25.0
25.0
25.0
25.0
25.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
74.0
75.0
75.0
75.0
75.0
75.0
75.0
75.0
75.0
75.0
99.0
99.0
99.0
99.0
99.0
99.0
99.0
99.0
99.0
99.0
99.0
99.0
99.0
99.0
99.0
99.0
99.0
99.0
99.0
99.0
99.0
99.0
99.0
99.0
99.0
99.0
99.0
99.0
99.0
99.0
99.0
99.0
99.0
99.0
99.0
99.0
99.0
99.0
99.0
99.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
33.0
33.0
33.0
33.0
33.0
33.0
33.0
33.0
33.0
33.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
67.0
67.0
67.0
67.0
67.0
67.0
67.0
67.0
67.0
67.0
99.0
98.0
98.0
98.0
98.0
98.0
98.0
98.0
98.0
98.0
98.0
98.0
98.0
98.0
98.0
98.0
98.0
98.0
98.0
98.0
98.0
98.0
98.0
98.0
98.0
98.0
98.0
98.0
98.0
98.0
98.0
98.0
98.0
98.0
98.0
98.0
98.0
98.0
98.0
98.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
2.0
2.0
2.0
2.0
2.0
2.0
2.0
2.0
2.0
2.0
2.0
2.0
2.0
2.0
2.0
2.0
2.0
2.0
2.0
2.0
2.0
2.0
2.0
2.0
2.0
2.0
2.0
2.0
2.0
2.0
2.0
2.0
2.0
2.0
2.0
2.0
2.0
2.0
2.0
26.0
26.0
26.0
26.0
26.0
26.0
26.0
26.0
26.0
26.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
73.0
73.0
73.0
73.0
73.0
73.0
73.0
73.0
73.0
73.0
99.0
99.0
99.0
99.0
99.0
99.0
99.0
99.0
99.0
99.0
99.0
99.0
99.0
99.0
99.0
99.0
99.0
99.0
99.0
99.0
99.0
99.0
99.0
99.0
99.0
99.0
99.0
99.0
99.0
99.0
99.0
99.0
99.0
99.0
99.0
99.0
99.0
99.0
99.0
99.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
25.0
25.0
25.0
24.0
24.0
24.0
24.0
24.0
24.0
24.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
75.0
75.0
75.0
76.0
76.0
76.0
76.0
76.0
76.0
76.0
99.0
99.0
99.0
99.0
99.0
99.0
99.0
99.0
99.0
99.0
99.0
98.0
98.0
98.0
98.0
98.0
98.0
98.0
98.0
98.0
98.0
98.0
98.0
98.0
98.0
98.0
98.0
98.0
98.0
98.0
98.0
98.0
98.0
98.0
97.0
97.0
96.0
96.0
96.0
96.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
2.0
2.0
2.0
2.0
2.0
2.0
2.0
2.0
2.0
2.0
2.0
2.0
2.0
2.0
2.0
2.0
2.0
2.0
2.0
2.0
2.0
2.0
2.0
3.0
3.0
4.0
4.0
4.0
4.0
27.0
27.0
27.0
27.0
27.0
27.0
27.0
27.0
27.0
27.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
73.0
73.0
73.0
73.0
73.0
73.0
73.0
73.0
73.0
73.0
99.0
99.0
99.0
99.0
99.0
99.0
99.0
99.0
99.0
99.0
99.0
99.0
98.0
98.0
98.0
98.0
98.0
98.0
98.0
98.0
98.0
98.0
98.0
98.0
98.0
98.0
98.0
98.0
98.0
98.0
98.0
98.0
98.0
98.0
98.0
98.0
98.0
98.0
98.0
98.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
2.0
2.0
2.0
2.0
2.0
2.0
2.0
2.0
2.0
2.0
2.0
2.0
2.0
2.0
2.0
2.0
2.0
2.0
2.0
2.0
2.0
2.0
2.0
2.0
2.0
2.0
2.0
2.0
34.0
34.0
34.0
33.0
33.0
33.0
33.0
33.0
32.0
31.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
66.0
66.0
66.0
67.0
67.0
67.0
67.0
67.0
68.0
69.0
simulations = repeat_simulations(100, 1000, InfectionRecovery(0.02, 0.002), 20)In the cell below, we plot the evolution of the number of alpha=0.5 inside the plot command).
let p = plot() for sim in simulations plot!(p, 1:1000, sim.I, alpha=.5, label=nothing) end # convert simulations (named tuple of arrays into matrix of infections) n_sims, n_time = length(simulations), length(simulations[1].I) i_matrix = zeros((n_sims, n_time)) for s in 1:n_sims i_matrix[s, :] = simulations[s].I # infection matrix end # cacluate the mean for each time period mean_infections = [mean(i_matrix[:,t]) for t in 1:n_time] plot!(p, 1:1000, mean_infections, lw=3, label="mean infections") pendxxxxxxxxxxplot_repeated_sims(simulations)plot_repeated_sims (generic function with 1 method)xxxxxxxxxxfunction plot_repeated_sims(simulation) n_time, n_sims = size(first(simulation).S)..., size(simulation)... p = plot() for sim in simulation plot!(p, 1:n_time, sim.I, alpha=0.5, label=nothing) end i_mean = mean(reduce(hcat,[sim.I for sim in simulation]), dims=2) plot!(p, i_mean, lw=3, label="average infections")end👉 Calculate the mean number of infectious agents of our simulations for each time step. Add it to the plot using a heavier line (lw=3 for "linewidth") by modifying the cell above.
Check the answer yourself: does your curve follow the average trend?
Hint
This exercise requires some creative juggling with arrays, anonymous functions, maps, or whatever you see fit!
xxxxxxxxxxmd"""👉 Calculate the **mean number of infectious agents** of our simulations for each time step. Add it to the plot using a heavier line (`lw=3` for "linewidth") by modifying the cell above. Check the answer yourself: does your curve follow the average trend?$(hint(md"This exercise requires some creative juggling with arrays, anonymous functions, `map`s, or whatever you see fit!"))"""👉 Write a function sir_mean_plot that returns a plot of the means of
xxxxxxxxxxmd"""👉 Write a function `sir_mean_plot` that returns a plot of the means of $S$, $I$ and $R$ as a function of time on a single graph."""sir_mean_plot (generic function with 1 method)x
function sir_mean_plot(simulations::Vector{<:NamedTuple}) n_sims, n_time = length(simulations), length(first(simulations).S) # convert simulations (named tuple of arrays into matrix of infections) s_mtx, i_mtx, r_mtx = zeros((n_sims, n_time)), zeros((n_sims, n_time)), zeros((n_sims, n_time)) for s in 1:n_sims s_mtx[s, :] = simulations[s].S # susceptible matrix i_mtx[s, :] = simulations[s].I # infection matrix r_mtx[s, :] = simulations[s].R # reinfection matrix end # calculate the mean for each time period s_mean, i_mean, r_mean = zeros(n_time), zeros(n_time), zeros(n_time) for t in 1:n_time s_mean[t], i_mean[t], r_mean[t] = mean(s_mtx[:,t]), mean(i_mtx[:,t]), mean(r_mtx[:,t]) end p = plot() plot!(p, 1:1000, s_mean, alpha=0.5, label="mean susceptible") plot!(p, 1:1000, i_mean, alpha=0.5, label="mean infected") plot!(p, 1:1000, r_mean, alpha=0.5, label="mean recovered") return pendxxxxxxxxxxsir_mean_plot(simulations)👉 Allow
xxxxxxxxxxmd"""👉 Allow $p_\text{infection}$ and $p_\text{recovery}$ to be changed interactively and find parameter values for which you observe an epidemic outbreak."""infection rate
recovery rate
x
begin simulations2 = repeat_simulations(1000000, 1000, InfectionRecovery(p_infection, p_recovery), 20) sir_mean_plot(simulations2)end👉 Write a function sir_mean_error_plot that does the same as sir_mean_plot, which also computes the standard deviation yerr=σ in the plot command; use transparency.
This should confirm that the distribution of
sir_mean_error_plot (generic function with 1 method)x
function sir_mean_error_plot(simulations::Vector{<:NamedTuple}) n_sims, n_time = length(simulations), length(first(simulations).S) # convert simulations (named tuple of arrays into matrix of infections) s_mtx, i_mtx, r_mtx = zeros((n_sims, n_time)), zeros((n_sims, n_time)), zeros((n_sims, n_time)) for s in 1:n_sims s_mtx[s, :] = simulations[s].S # susceptible matrix i_mtx[s, :] = simulations[s].I # infection matrix r_mtx[s, :] = simulations[s].R # reinfection matrix end # calculate the mean for each time period s_mean, i_mean, r_mean = zeros(n_time), zeros(n_time), zeros(n_time) s_stdev, i_stdev, r_stdev = zeros(n_time), zeros(n_time), zeros(n_time) for t in 1:n_time s_mean[t], i_mean[t], r_mean[t] = mean(s_mtx[:,t]), mean(i_mtx[:,t]), mean(r_mtx[:,t]) s_stdev[t], i_stdev[t], r_stdev[t] = Statistics.std(s_mtx[:,t]), Statistics.std(i_mtx[:,t]), Statistics.std(r_mtx[:,t]) end p = plot() plot!(p, 1:1000, s_mean, alpha=0.25, ribbon=s_stdev, label="mean S") plot!(p, 1:1000, i_mean, alpha=0.25, ribbon=i_stdev, label="mean I") plot!(p, 1:1000, r_mean, alpha=0.25, ribbon=r_stdev, label="mean R") return pendsir_mean_error_plot2 (generic function with 1 method)xxxxxxxxxxfunction sir_mean_error_plot2(simulations::Vector{<:NamedTuple}) n_sims = length(simulations) xs = 1:1000 # mtx_s is size(1000, 20) - rows are time, cols are simulations mtx_s = reduce(hcat,[sim.S for sim in simulations]) mtx_i = reduce(hcat,[sim.I for sim in simulations]) mtx_r = reduce(hcat,[sim.R for sim in simulations]) #sum across dims=2 collapses the 2nd dim resulting in size(1000,1) array mean_s = sum(mtx_s, dims=2) ./n_sims mean_i = sum(mtx_i, dims=2) ./n_sims mean_r = sum(mtx_r, dims=2) ./n_sims std_s = (sum((mtx_s .-mean_s).^2, dims=2) ./ (n_sims-1)) .^0.5 std_i = (sum((mtx_i .-mean_i).^2, dims=2) ./ (n_sims-1)) .^0.5 std_r = (sum((mtx_r .-mean_r).^2, dims=2) ./ (n_sims-1)) .^0.5 p = plot() plot!(p, xs, mean_s, alpha=0.25, ribbons=std_s, label="mean S") plot!(p, xs, mean_i, alpha=0.25, ribbons=std_i, label="mean I") plot!(p, xs, mean_r, alpha=0.25, ribbons=std_r, label="mean R") return pendsir_mean_error_plot2(simulations2)xxxxxxxxxxsir_mean_error_plot(simulations2)Exercise 3.3
👉 Plot the probability distribution of num_infected. Does it have a recognisable shape? (Feel free to increase the number of agents in order to get better statistics.)
let N = 100 # number of agents T = 1000 # number of sweeps agents = generate_agents(N) infection = InfectionRecovery(p_infection, p_recovery) for i ∈ 1:T sweep!(agents, infection) end num_infected = [agent.num_infected for agent in agents] p = bar(frequencies(num_infected)) vline!(p, [mean(num_infected)], label="mean")endExercse 3.4
👉 What are three simple ways in which you could characterise the magnitude (size) of the epidemic outbreak? Find approximate values of these quantities for one of the runs of your simulation.
xxxxxxxxxxmd"""1. Average number of people infected (how contagious) -> 0.952. Peak infection rate and expected time of peak => 65% @ T = 9003. Peak death rate and expected time of peak"""Exercise 4: Reinfection
In this exercise we will re-use our simulation infrastructure to study the dynamics of a different type of infection: there is no immunity, and hence no "recovery" rather, susceptible individuals may now be re-infected
Exercise 4.1
👉 Make a new infection type Reinfection. This has the same two fields as InfectionRecovery (p_infection and p_recovery). However, "recovery" now means "becomes susceptible again", instead of "moves to the R class.
This new type Reinfection should also be a subtype of AbstractInfection. This allows us to reuse our previous functions, which are defined for the abstract supertype.
xxxxxxxxxxmutable struct Reinfection <:AbstractInfection p_infection p_recoveryend👉 Make a new method for the interact! function that accepts the new infection type as argument, reusing as much functionality as possible from the previous version.
Write it in the same cell as our previous interact! method, and use a begin block to group the two definitions together.
Exercise 4.2
👉 Run the simulation 20 times and plot
Note that you should be able to re-use the sweep! and simulation functions , since those should be sufficiently generic to work with the new step! function! (Modify them if they are not.)
x
let A = 100 # number of agents N = 20 # number of trials T = 1000 # re_sim = repeat_simulations(A, T, Reinfection(0.02, 0.002), N) plot_repeated_sims(re_sim) #md""" infections plateu but do not decline"""end👉 Run the new simulation and draw
x
let A = 100 # number of agents N = 20 # number of trials T = 1000 # re_sim = repeat_simulations(A, T, Reinfection(0.02, 0.002), N) sir_mean_error_plot(re_sim) #md""" infections plateu but do not decline"""endExercise 5 - Lecture transcript
(MIT students only) Please see the link for hw 4 transcript document on Canvas. We want each of you to correct about 400 lines, but don’t spend more than 15 minutes on it. See the the beginning of the document for more instructions. :point_right: Please mention the name of the video(s) and the line ranges you edited:
Abstraction, lines 1-219
Array Basics, lines 1-137
Course Intro, lines 1-44
(for example)
xxxxxxxxxxlines_i_edited = md"""Abstraction, lines 1-219Array Basics, lines 1-137Course Intro, lines 1-44 (_for example_)"""Before you submit
Remember to fill in your name and Kerberos ID at the top of this notebook.
Function library
Just some helper functions used in the notebook.
hint (generic function with 1 method)almost (generic function with 1 method)still_missing (generic function with 2 methods)keep_working (generic function with 2 methods)Fantastic!
Splendid!
Great!
Yay ❤
Great! 🎉
Well done!
Keep it up!
Good job!
Awesome!
You got the right answer!
Let's move on to the next section.
correct (generic function with 2 methods)not_defined (generic function with 1 method)